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Abstract. - Porous sediments in geological systems are exposed to stress by the above-laying mass 
and consequent compaction, which may be significantly nonuniform across the massif. We derive 
scaling laws for the compaction of sediments of similar geological origin. With these laws, we eval- 
uate the dependence of the transport properties of a fluid-saturated porous medium (permeability, 
effective molecular diffusivity, hydrodynamic dispersion, electrical and thermal conductivities) on 
its porosity. In particular, we demonstrate that the assumption of a uniform geothermal gradient 
is not adequate for systems with nonuniform compaction and show the importance of the derived 
scaling laws for mathematical modelling of methane hydrate deposits; these deposits are believed 
to have potential for impact on global climate change and Glacial-Interglacial cycles. 
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Introduction. — The reconstruction of properties of 
grounds is an important problem related to geological sur- 
veys for extraction of minerals and hydrocarbons, con- 
struction of buildings, forecasting geological hazards (seis- 
mic, erosion-related, anomalous impact on the climate, 
etc.) [TH3]. Dealing with this problem one faces challenges, 
some of which hardly may be overcome. These challenges 
are related to the impossibility of direct measurements 
of required parameters across large massifs. Even mak- 
ing boreholes provides limited information about the nar- 
row vicinity of the borehole; for instance, measurements of 
the electrical conductivity and the porosity of the porous 
medium are generally not enough to reconstruct its per- 
meability, which practically can not be measured directly. 
Thus, the problem actually turns into one of the recov- 
ery of relations between different (transport) parameters 
of the porous media. Generally, this problem is non- 
resolvable, because it requires thorough knowledge of the 
composition of the massif, its geological and seismic his- 
tory, etc. Meanwhile, many recent studies deal with sys- 
tems where the massif possesses a homogenous geologi- 
cal origin on the field scale [2HH]- Opportunities for an 
advance in the problem of reconstruction of relations be- 
tween parameters for such kinds of systems might lay in 
the field of mathematical physics. In this study we wish to 
approach this problem in application to some important 



geological systems, like ocean bed with methane hydrate 
deposits. 

The ocean bed in regions with intense mud rain is very 
attractive and important for research due to bio- and geo- 
logical richness and activity [4][5] ; for example, these ocean 
bed systems host marine methane hydrate deposits. In 
such an ocean bed, sediments are exposed to a pressure 
load and have experienced a certain history of this load. 
These two factors result in compaction of the porous sedi- 
ments. A typical sample of such a compaction, increasing 
with depth, can be seen in Fig.[T^, where the sediment 
porosity ((/>) for different depths below the water-sediment 
interface is reported [H[S]. According to [51[7] one can fit 
the observed porosity with the exponential function 



(z) = o exp(-z/L), 



(1) 



where z is the depth below the water-sediment interface, 
L is the characteristic depth of compaction (see Caption 
to Fig.CTJ). 

Earlier approaches to the relation between porosity and 
transport properties (such as, e.g., the Kozeny-Carman 
relation for flow [31IID]) adopted additional (though rea- 
sonable) constrains on the geometry of pore channels 
and utilized several characteristic parameters (for review 
see [3]). As well as the family of empirical equations re- 
lating permeability to porosity and irreversible water sat- 
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Fig. 1: (a) Observed porosity vs the depth below the water- 
sediment interface for the ODP (Ocean Drilling Program 4,5) 
site 997 on the Blake Ridge crest — one of the largest marine 
hydrate provinces. The solid line represents the exponential de- 
crease of the porosity (Eq. |l} with 4>o = 0.69 and L — 2km [6]). 
(b) Predicted permeability [blue dashed line, Eq. Q] and ther- 
mal conductivity [red dash-dot line, Eq. (|10|l ] of sediments. 



uration these models do not provide information on 
how the characteristic parameters will be transformed un- 
der compaction. Thus, a particular theory of the latter is 
remaining required. 

In this Letter, in order to reconstruct the relations be- 
tween parameters, we first argue for geometric features of 
compaction and corresponding scaling laws. Then we de- 
rive the dependence of transport coefficients on the poros- 
ity and, for a particular example, apply the derived de- 
pendence to examine whether the assumption of a linear 
temperature profile in sediments is adequate. We find that 
the amended predicted rate of the production of methane 
hydrate differs significantly from the "traditional" predic- 
tions assuming the linear temperature profile; this finding 
is of importance for the assessment of the planetary inven- 
tory of methane hydrate. 

Physical model and scaling laws: Compaction 
of porous medium. — Prior to constructing a mathe- 
matical description of the problem, we need to establish 
physical features of the sediment compaction in Fig.QJi. 

In general, when considering the variation of transport 
coefficients with porosity, one should keep in mind two 
critical thresholds for porosity [3]. First, there is a crit- 
ical porosity important for acoustic processes in porous 
materials; a more porous material can be treated as a 
"suspension" and a less porous one is "solid". This first 
critical porosity is material-specific and is typically near 
0.5, i.e., might be close to the lower edge of the part of the 
sediment column in Fig. [I] For porosities lower than the 
second threshold, chemical compaction is taking over the 



mechanical one. This second transition occurs for porosi- 
ties which are beyond the scope of our study. According 
to seismic data, for the part of the sediment column we 
consider, the first critical transition does not occur as well. 
Indeed, geological and mechanical structure of sediments 
is reported to be continuous enough for using the seismic 
wave reflection to detect as small volumetric fraction of 
bubbles in pores as 1-2% [5l[T2j[T3]. Hence, qualitative 
transitions in the structure of porous sediments are not 
expected. 

The stress load on the system is anisotropic; the ver- 
tical direction is discriminated due to gravity. In a solid 
body the stress would be strictly anisotropic. However, 
the stress in granular materials (including cemented ones) 
is known to be not distributed homogeneously but to 
form force chains [13]. The stress along chains is much 
bigger than the average stress, and the grain displace- 
ments and deformations in the course of compaction are 
mainly driven by these chains. Meanwhile, the random- 
ness of the geometry of contacts and the branching of 
force chains decreases the anisotropy of the network of 
force chains in comparison to the stress anisotropy for a 
solid body [TH - fll)] . For simplicity, we assume the com- 
paction process to be isotropic. Henceforth, we treat an 
isotropic compaction and consider the sediments to have 
similar geologic origin (which is reasonable for the geologi- 
cal systems we consider, on the timescale of several million 
years [5]). 

In our model, as a first approach to the problem, we 
do not consider the fragmentation of sediment grains. 
Anisotropy of compaction and fragmentation are the main 
reasons for the change of the topology of the pore channels' 
network. Without them, compaction affects this topology 
non-efficiently, but rather shrinks the channels. Further, 
the variation of the density of the solid matrix material is 
negligible against the background of the pore shrinking. 

Let us derive general scaling laws for the compaction 
of sediments the physical features of which are described 
above. We consider a certain volume V of the porous 
medium and fix this volume to particles of the solid ma- 
trix (V changes due to compaction); L, I, and V s stand for 
the total channels' length in V, the characteristic channels' 
diameter, and the volume of the solid material in V, re- 
spectively. Hence, 
(i) porosity 



V 



LI 2 



V s + Ll' 2 



and thus 



(ii) Owing to the unchanging shape of the channels' net- 
work, L 3 oc V . 

(iii) Owing to the unchanging density of the solid mate- 
rial, V s = (1 — <j))V is constant in the course of compaction; 
therefore V (x (\ — <fi)~ 1 and 



L oc 



(1-^)V3- 



(2) 
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From (i), 



I i 

~>L ^ (1-0)2/3- 



(3) 



These relations are obviously relevant for a tube model 
of the porous medium. However, they are derived without 
actual relaying on features of any specific model and we 
expect them to be reasonably accurate for realistic geo- 
metric models of porous medium [17j as long as we deal 
with a sparse porous structure (in Fig. Q] porosity varies in 
the range from 0.5-0.7). The scaling of the linear measure 
I for realistic models describes the variation of transversal 
linear measures of pores. For low-porosity materials the 
actual scaling laws become sensitive to features of the pore 
geometry (e.g., see [17]). 

Permeability. We use the definition of the permeabil- 
ity K according to the Darcy's law (e.g., see ffilTB"]): 

v f = -(K/rj)Vp, 

where Vf is the filtration velocity of the fluid, r\ is the 
dynamic viscosity, p is pressure. When the distribution 
of the orientation of pore channels and topology of the 
pore network are unchanged — as they are in our model of 
compaction — the fluid speed in pores Uf = 4>~ 1 Vf requires 
the pressure gradient proportional to Z~ 2 . Thus, Vf oc 
(j)l 2 \Vp\ and 



K oc 



(1 -0) 2 /3- 



(4) 



For instance, the data reported in Fig.[T] yield K(<f> = 
O.69)/K(0 = 0.49) w 2.8 due to compaction; that is the 
permeability of the upper sediment zone (0 « 0.69) is by 
factor 2.8 larger than that at the bottom of the shown 
sediment column (0 w 0.49). 

Remarkably, in [§] a rough dependence of the perme- 
ability on the porosity, K oc 2 [cf Eq. ((3])], was adopted 
because of the lack of reasonable theories on scaling laws 
for sediments experiencing compaction. 

Molecular diffusion of solute. The evolution of the so- 
lution concentration C in quiescent pore water is governed 
by the equation 



^(0C) = V-[ 7 z 5 (0)^ m VC], 



(5) 



where D m is the molecular diffusivity of the solute in bulk, 
and 7d(0) is the geometric factor featuring the pore net- 
work. 

Similarly to the permeability, the geometric factor for 
the molecular diffusivity (jd) depends on porosity which 
varies significantly with depth. The length of the channels 
L does not effect the diffusional flux until the statistics of 
channel orientations are changed. With the concentra- 
tion gradient given, the solute flux through the area S is 
linearly proportional to the area of channel cross-section, 

Spore- 

Id oc S porc /S = 0. (6) 



Recall, porosity is exactly the average value of S polc /S. 
This becomes evident if one considers the cubic volume 
thinly sliced parallel to one of its sides; for each slice 
the fraction of the pore volume is Spore/ S, and, thus, for 
the whole cube volume the ratio of pore volume to the 
cube volume, which is porosity 0, is the average value of 

Spore/ S . 

Hydrodynamic dispersion. Due to the irregularity of 
the microstructure of the pore network, the macroscopi- 
cally uniform displacement of liquid results in a mixing 
flow in pores, which acts as an additional diffusion and is 
referred to as hydrodynamic dispersion [TJ[T5]. The hy- 
drodynamic dispersion in an isotropic medium is strictly 
anisotropic; the longitudinal and lateral dispersion coeffi- 
cients differ and are linearly proportional to the filtration 
speed Vf [15] : 

D|| = Vfdi, D± = vfd 2 . 

For a steady viscous flow in pores, D^± oc u 2 T corr oc 
Uflcorr (recall, the fluid speed in pores Uf = Vf /4>). For 
an isotropic compaction, i corr is scaled as the pore network 
skeleton, that is oc L, and Eq. §3§ yields -Dim oc (vf /<j))L oc 
v f /[<j)(l - 0) 1/3 ]. Hence, 



1 



0(1-0)V3- 



(7) 



Notice, the geometric factor j£> for molecular diffusion 
[Eq. ©] and the hydrodynamic dispersion coefficients 
[Eq. ([7])] are affected by the compaction differently. 

Electrical conductivity. For the electrical conductiv- 
ity one should clearly distinguish two cases: (i) sea water 
and (ii) pure water in pores. Due to electrolytes dissolved, 
sea water possesses the electrical conductivity about 5S/m 
which is much more than that of the porous-skeleton ma- 
terial. The current flows through the liquid volume. On 
the pore scale, this case is geometrically equivalent to the 
case of molecular diffusion. Indeed, for a steady diffusive 
flux we have the equation AC = for concentration C 
in the pore volume, zero normal derivative of C on the 
pore walls, and fixed mean (macroscopic) gradient of C; 
for electrostatic potential if we find the same equation 
Aip = in the pore volume with zero normal derivative of 
ip on the pore walls and fixed macroscopic gradient of ip. 
For the electrical conductivity a of sea water, this equiv- 
alence yields 

a oc , (8) 

the same scaling law as Eq. ^ . This law can be observed 
f° r 4> si 0.15 in realistic models of the pore morphology 
and experiments (see Ref. [17] and references therein for 
experimental data). 

The case of pure water is significantly more subtle. 
Without electrolytes dissociated, water has small num- 
ber of charge carriers; the mineral surface conductivity 
can make significant contribution into the microscopic 
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electrical conductivity. For sands, experiments demon- 
strate nearly the same conductivity for wet massif and 
the massif fully saturated with pure water, which indi- 
cates that the electrical current flows along the water- 
mineral interface [HI. In this case, resistivity is mainly 
contributed by the sharp sand grain contacts; the geom- 
etry of these contacts is controlled by many factors, in- 
cluding the stress [21] . This problem is very complex and 
lies beyond the immediate applicability field of our results. 
Thus, we emphasize that Eq. (jHJ is valid for salt water only. 

Heat diffusion. Heat transfer in the porous medium 
is governed by the equation 



under time-independent conditions. Hence, 



4>)p s cp.s + 4>PfCp,f)T] 
+V • [vfPfCpjT] = V ■ K0)VT] , 



(9) 



where p s , pf and cp. s , cpj are the densities and the spe- 
cific heat capacities of the solid matrix and the fluid, 
respectively, k{4>) is the heat conductivity of the fluid- 
saturated medium. 

An evaluation of the dependence of the macroscopic 
thermal conductivity on the porosity of sediments expe- 
riencing compaction is much more complicated than for 
the permeability and the solution diffusivity, because heat 
flows through both the solid matrix and the fluid in pores 
and the fluxes in these two subsystems (with complex ran- 
dom geometry) should be found and conjuncted at the in- 
terface. This issue requires a particularly accurate study. 
Straightforward scaling rules cannot be derived with our 
approach here and, instead, we relay on empirical relation 
suggested in Ref. 22 for porous media saturated with a 
weakly heat conducting fluid (such as air of water, which 
has 5-10 times smaller heat conductivity than typical min- 
eral materials). 

Chaudhary and Bhandari [32] reported the law 



k(4>) = k s — 



4> + (1 - 4>)Kf/K s 



(10) 



where Kf and the conductivities of matrix mate- 

rial and fluid in pores, respectively, and n — 0.5(1 — 
In 4>) / In[</>(1 — 0)k s /k/], to be satisfactory accurate for 
a broad variety of sediment-kind porous materials. For 
our study of carbon-bearing sediments, the relevant ma- 
trix material heat conductivity is k s — 2.93 J/(msK) [22., 
(for water Kf = 0.58J/(msK) is well known). In partic- 
ular, for the porosity profile in Fig. [I] n{<j) — 0.69)/k(</> = 
0.49) ~ 0.66, i.e., the temperature profile slope varies by 
factor 1.5 for different parts of the sediment column. 

Geothermal gradient and methane hydrate in- 
ventory. — While the researchers modelling methane 
hydrate deposits {e.g., [BHS]) adopt a uniform geothermal 
gradient G := dT/dz, it is, in fact, nonuniform. Instead, 
the heat flux, which is the product of the geothermal gradi- 
ent and the (nonuniform) thermal conductivity, is uniform 



G(z) 



dT{z) [heat flux] k{0) dT(0) 



dz 



k(z) 



k(z) dz 



Eq. CETJ yields 



T(z) = T(0) + Go 



«(0) 
k(zi) 



dzi, 



(11) 



(12) 



where Go is the geothermal gradient next to the water- 
sediment interface. The latter expression is convenient 
when Go is directly measured in shallow upper layer of 
sediments. However, for sediments bearing methane hy- 
drate, the geothermal gradient is practically derived in a 
different way. 

Immediately below the floor of deep seas the pressure of 
the water column is large enough and the temperature is 
low enough for methane hydrate to be stable. Deeper into 
sediments the temperature grows and, at certain depth, 
the pressure becomes not sufficiently large to stabilize 
the hydrate (the critical pressure depends on temperature 
nearly exponentially [23]). Thus, the Hydrate Stability 
Zone (HSZ) spreads in sediments from the water-sediment 
interface down to a certain depth, the base of the HSZ, 
which is detected by the reflection of seismic waves [12] . 
Practically, for marine methane hydrates the geothermal 
gradient is inferred from the position z = Lh of the base 
of the HSZ [H]: G := (T(L h )-T(Q))/L h , where T{L h ) is 
calculated as the hydrate destabilization temperature for 
the known water salinity and hydrostatic pressure P{Lh) 
of the water column {e.g., one can employ an accurate 
mathematical model of hydrate from [24). For this case, 

T{z) = T(0) + (T{L h ) - T(0)) £*~y • (13) 

J k L {zi)dzi 

In Fig. [5^, the temperature profile consistent with com- 
paction is compared to "traditional" linear profiles guessed 
either from the observed position of the base of HSZ or 
from measurements of the sediments temperature profile 
next to the water-sediment interface. Remarkably, assum- 
ing the geothermal gradient being uniform is especially 
inaccurate for the latter case (red dotted line): the as- 
sumption of the linear profile significantly rises the base 
of the HSZ. 

The role of the nonuniformity of G appears to be es- 
pecially pronounced in the problem of hydrate forma- 
tion. When a hydrate is present in the HSZ, (i) the 
concentration of methane dissolved in water is deter- 
mined by the thermodynamic equilibrium between the hy- 
drate and the aqueous solution, in other words, it equals 
the solubility, and (ii) the diffusive flux of the aqueous 
methane may posses non-zero divergence equal to the for- 
mation/dissociation rate of hydrate in pore water (up to 
a known constant multiplier determining the fraction of 
methane in hydrate). We employed the thermodynamic 
model of a hydrate developed in [21] for the calculation of 
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Fig. 2: (Color on-line) Predictions for the data in Fig.[T] Temperature profiles (a), the aqueous methane solubility in equilibrium 
with hydrate (b), and the consequent production of methane hydrate due to diffusion of methane dissolved in pore water (c) 
for an accurate nonuniform geothermal gradient (black solid lines), and for simplified linear temperature profiles guessed either 
from the observed position of the base of HSZ (blue dashed lines) or from the observed geothermal gradient immediately next 
to the water-sediment interface (red dotted lines). 



the solubility profiles (Fig.[5jD) and then derived the con- 
tribution of this flux into the hydrate production (Fig. [2}:) . 
In addition, to calculate the divergence of the diffu- 
sive flux we accounted for the strong dependence of the 
methane molecular diffusivity on the temperature (the de- 
pendence on pressure is negligible): D m w (7+0.4K~ 1 (T— 
273. 15K) + 1CT 3 K- 2 (T - 273. 15K) 2 ) • ICT 10 m 2 /s, which 
fits well with experimental data {e.g., |25j). One can see, 
that linear temperature profiles significantly overestimate 
the production of methane hydrate in the lower part of the 
HSZ. Thus, the results of mathematical modelling which 
simultaneously considers the compaction of the sediment 
and ignores the consequent non-linearity of the tempera- 
ture profile {e.g., [EH9]) are significantly affected by this 
inconsistency. 

Notice, here we argue for the importance of the scal- 
ing laws for compaction and assess the physical accuracy 
of models adopted in [BH9] in this context only. This is 
the reason, why we readdress the pure Fickian diffusion of 
aqueous methane and do not consider non-Fickian effects, 
thermodiffusion and gravitational stratification of the so- 
lute, ignored in the existing models of the formation of 
hydrate deposits, although the importance of non-Fickian 
effects was shown in [2"rJII27| . 

Conclusion. — Summarizing, we have described an 
isotropic compaction of porous medium and derived scal- 
ing laws for geometrical properties of the pore structure. 
These laws have yielded dependencies of transport prop- 
erties [permeability, Eq. Q , effective molecular diffusivity, 
Eq. dHl), hydrodynamic dispersion, Eq. (J7J), electrical con- 
ductivity, Eq. ((U, and thermal conductivity, Eq. (ff"0|) ] on 
the porosity for porous sediments of similar geological ori- 
gin. Notably, the compaction of sediments (for example, 
see Fig.QJi) is an inherent feature of most geological sys- 
tems on the field scale. In particular, for paradigmatic 
models of formation of marine methane hydrate [6 9 , 



compaction is a "key ingredient" . The employment of 
our results for transport coefficients provides an oppor- 
tunity for a significant enhancement of physical soundness 
and relevance of the modelling of sediments experiencing 
compaction and, in particular, the global methane hydrate 
inventory {e.g., Fig.[5J; demonstrates the inaccuracy of the 
hydrate production rate in treatments disregarding the 
variation of the thermal conductivity due to compaction) . 
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